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[i]  We  consider  diffusion  of  outer  zone  radiation  belt 
electrons  by  chorus  waves.  Quasi-linear  diffusion 
coefficients  valid  outside  the  plasmasphere  have  only  been 
calculated  recently,  and  indicate  that  the  energy  and  cross 
diffusion  rates  can  be  comparable  to  that  for  pitch  angle 
diffusion.  Proper  solution  of  the  diffiision  equation  for 
phase  space  density  must  therefore  be  based  on  the  full 
diffusion  tensor,  but  this  has  been  plagued  by  numerical 
problems  associated  with  the  large  and  rapidly  varying 
cross  terms.  To  circumvent  this,  techniques  are  developed 
for  transforming  to  variables  in  which  the  cross  diffusion 
term  vanishes.  A  model  calculation  shows  significant 
diffusion  of  phase  space  density  at  Z,  =  4.5  from  0.2  MeV 
up  to  a  few  MeV  in  less  than  a  day.  Citation:  Albert,  J.  M., 
and  S.  L.  Young  (2005),  Multidimensional  quasi-linear  diffusion 
of  radiation  belt  electrons,  Geophys.  Res.  Lett.,  32,  L14110, 
doi;  1 0. 1 029/2005GL023 191. 


1.  Introduction 

[2]  A  standard  approach  to  studying  the  dynamics  of 
energetic  particles  in  Earth’s  radiation  belts  is  via  a  quasi- 
linear  diffosion  equation,  which  describes  the  phase  space 
distribution  function  suitably  averaged  over  gyro-,  bounce, 
and  (longitudinal)  drift  frequencies.  Electromagnetic  waves 
resonant  with  these  frequencies  lead  to  diffusion  in  the  three 
adiabatic  invariants  associated  with  them;  other  terms  due 
to,  e.g..  Coulomb  collisions  are  easily  incorporated  as  well. 
In  principle  coupling  occurs  between  the  diffusion  in  all 
three  invariants,  but  in  a  nearly  axisymmetric  dipole  the 
diffusion  in  the  third  invariant  J3  (associated  with  waves  at 
the  drift  frequency)  is  largely  decoupled  from  diffusion  in  Ji 
and  J2  (driven  by  waves  resonant  near  or  below  the  much 
higher  cyclotron  frequency).  However,  the  diffusion  in 
and  J2  is  highly  coupled,  according  to  the  formulation  of 
Lyons  [1974a,  1974b].  After  bounce  averaging,  this  treat¬ 
ment  gives  diffusion  coefficients  for  equatorial  pitch  angle 
oo,  momentum  (or  energy),  and  cross  diffusion  D^^p,  which 
can  be  recast  in  terms  of  the  adiabatic  invariants  [Schulz  and 
Lanzerotti,  1974].  Then  the  diffusion  equation  takes  the 
form 


dt 


d 

dJi 


(r.  ^ 

Dn  p7-  +  Z3i2 
^  oJ\ 


+ 


dJi 


(A2g  +  Z)22^) 


9  ^  df 


(1) 


[3]  Inside  the  plasmasphere,  waves  cyclotron  resonant 
with  radiation  belt  electrons  occur  primarily  as  whistler 
mode  hiss.  There  the  plasma  frequency  fpe  is  large  compared 
to  the  cyclotron  frequency  fee,  resulting  in  a  pitch  angle 
diffusion  rate  much  greater  than  the  energy  diffusion  rate, 
with  the  cross  diffusion  rate  intermediate  between  the  two 
[Lyons,  1974b].  Thus  pure  pitch  angle  diffusion  combined 
with  radial  diffusion  is  a  reasonable  approximation  there; 
further  simplifying  the  pitch  angle  diffusion  to  a  loss  rate 
yields  a  successful  description  of  the  inner  electron  radiation 
belts,  at  least  during  quiet  or  moderately  disturbed  condi¬ 
tions  [Lyons  and  Thome,  1973].  Outside  the  plasmasphere, 
however,  such  models  have  not  been  able  to  reproduce  the 
localized  high  energy  phase  space  density  peaks  frequently 
observed  to  develop  during  the  recovery  phase  of  magnetic 
storms  [e.g.,  Brautigam  and  Albert,  2000;  Shprits  and 
Thome,  2004].  It  has  been  argued  that  the  condition  fpe  ~ 
fee  enables  storm-time  enhanced  whistler  mode  choms  to 
drive  rapid  electron  acceleration  [Summers  et  al,  1998; 
Home  etal,  2003, 2005;  Albert,  2005],  which  could  account 
for  the  energized  population.  Although  the  development  of 
chorus  waves  involves  strongly  nonlinear  processes  [Sazhin 
and  Hayakawa,  1992],  we  assume  that  quasi-linear  theory  is 
a  meaningful  description  of  the  effect  of  fully  developed 
chorus  on  the  high  energy  population.  Suggestive  results 
have  been  obtained  from  a  one  dimensional  energy  diffusion 
equation  averaged  appropriately  over  oo  [Summers  and  Ma, 
2000;  Summers  et  al.,  2004;  Home  et  al,  2005;  Li  et  al, 
2005]  but  studying  this  process  in  detail  requires  treating  the 
full  velocity  space  diffusion  equation.  However,  straightfor¬ 
ward  solution  of  equation  (1)  is  elusive,  even  when  the  usual 
conditions  for  local,  linear  numerical  stability  are  well 
satisfied,  because  of  the  large  and  rapidly  varying  cross 
diffusion  coefficient.  (See  Albert  [2004]  for  a  review  of 
previous  approaches.)  To  overcome  this,  new  variables  are 
constmeted  for  which  the  cross  diffusion  term  vanishes,  and 
standard  finite  difference  techniques  are  applied  to  the 
transformed  diffusion  equation. 


2.  Eliminating  the  Cross  Term 

[4]  In  general,  when  transforming  from  {Jx,  J2)  to  new 
variables  {Qi,  Q2),  the  diffusion  matrix 


Djj  = 


Al 

A2 


A2 

D22 


(2) 


becomes  [Haerendel,  1968;  Schulz,  1991] 


Because  of  the  decoupling,  the  (or  L  or  “radial”) 
diffusion  is  relatively  easy  to  treat,  even  if  it  is  large. 


dQxIdJx 

dQi/dJx 


dQxlQJi^r,  \9QildJx 
dQi/dJi  dQi  /dJ2  dQi/dJi 


(3) 
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It  is  desired  that  the  cross  diffusion  terms  in  the  new 
variables  vanish.  It  will  be  assumed  that  the  diffusion  is 
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[6]  A  second  condition  is  required  in  order  to  determine 
the  two  slopes.  One  alternative  is  to  impose  ^2  =  —  1,  so 
that  the  families  of  constant-gj  and  constant-g2  curves  are 
orthogonal  in  the  (Ji,  plane;  this  is  equivalent  to  taking 
the  constant-g  curves  parallel  to  the  eigenvectors  of  the  Z)jrj 
matrix.  A  simpler  approach,  adopted  here,  is  to  specify  gi 
as  a  convenient  function  of  Ji  and  J2,  which  determines 
and  thus  S2  from  equation  (5). 

[7]  As  a  simple  illustration,  take  the  elements  of  D jj  to  be 
constant.  The  Green  function  solution  of  the  (Jj,  J-^  part  of 
equation  (1),  for  a  distribution  initially  localized  at  (Jio,  /20) 
far  from  the  boundaries,  is 

exp[-(D22&/?  -  2D,26i/i&/2  +Dii&/|)/4tA]/4iU\/A, 


where  hJ j  —J\  —  ./io»  b/2  ~  J2  —  •^20*  and  A  —  D\ \D22  —  -^^2* 
Choosing  Qi  =  J2  leads  to  Q2-J\  —  (^>12/^22^2  and  G  = 
— 1,  as  well  as  D]  =  D22  and  D2  =  Z)ii  -  D12/D22.  Then 
it  may  be  checked  that  the  Green  fimction  solution  of 
equation  (7),  namely  exp(-6gi/4Z)if  -  hQ2lAD2t)l 
4'xr\/£)iZ)2,  is  the  same  as  the  expression  above. 

[s]  Choosing  Q\  to  be  the  equatorial  pitch  angle  0(o,  the 
equation  for  constant-g2  curves  can  be  written  as 


Figure  1.  First  three  panels:  inverse  timescales  from  quasi- 
linear  diffusion  coefficients,  in  units  of  s~’,  for  electrons  at 
L  =  4.5  based  on  a  global  model  of  storm  time  whistler¬ 
mode  chorus  waves.  The  last  panel  shows  where  the  cross 
diffusion  coefficient  is  positive  (red)  or  negative  (blue). 

truly  two  dimensional,  so  that  det(£)jy)  >  0.  This  is  strictly 
true  for  the  boimce-averaged  quasi-linear  diffusion  coeffi¬ 
cients  considered  here,  though  barely  so  if  the  resonances 
are  highly  localized  in  both  wave  normal  angle  and  latitude 
(as  can  occur  for  particles  with  «  =  0  resonances  only) 
[Albert,  2004]. 

[5]  The  slopes,  dJ2ldJ\,  of  curves  on  which  Qi  and  Q2 
are  constant  are  given,  respectively,  by 


dE  _  Pg^E 
dcto  Pctaao 

and  the  transformed  diffusion  coefficients  simplify  to 


The  (qo,  E)  diffusion  coefficients  and  {Du,  D22)  are 
related  by  equation  3  so  that,  e.g.,  Da^og  ~  (Aoior/At  rather 
than  (pAoioflAt  as  by  Lyons  [1974a,  1974b].  As  above,  the 
expressions  for  the  quasi-linear  diffusion  coefficients 
guarantee  that  £>2  is  positive. 


o  dQi/dA  dQ2/dJi 

'  dQi/dJi'  '  dQi/dJi 

From  equation  (3),  the  requirement  D12 
expressed  in  terms  of  the  slopes  as 

SiS2Du  —  (S'!  +  S2)Di2  +  D22  =  0,  (5) 

and  the  nonvanishing  diffiision  coefficients  are 

D,  =  {dQ,/dJ2Y  (S?D„  -25,D,2  -1-022), 

(6) 

O2  =  {dQ2/dJ2f  (5|0„  -2520,2  -I-O22), 

which  are  both  positive  since  O,,  D22  >  £>12  [Albert,  2004]. 
In  the  new  variables,  the  (J,,  J2)  part  of  the  diffusion 
equation  simplifies  to 


(4) 

=  0  can  be 


dt 


LfJ-rn  M. 

G  dQ\  dQ2  dQ2 


)■ 


(7) 


where  G  =  det  [d{J\,  J^ld{Q\,  and  there  is  no  cross 
diffusion  term.  The  radial  diffusion  term,  if  included,  is 
unchanged. 


3.  Performing  the  lYansformation 

[9]  The  diftoion  coefficients  {Dg^g^,  Dg^plp,  Dpplp^,  all 
with  units  of  s  ')  due  to  storm-time  chorus  were  calculated 
at  L  =  4.5,  using  the  wave  models  of  Home  et  al.  [2005]  and 
the  computational  methods  of  Albert  [2005].  In  Aese  wave 


0  15  30  45  60  75  90 

«o 


Figure  2.  Curves  of  constant  g,  =  oq  (dotted)  and  g2 
(solid)  corresponding  to  the  diffusion  coefficients  of 
Figure  1,  constmcted  so  that  the  cross  diffusion  coefficient 
vanishes. 
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□  i  D2 


0  30  60  90  0  30  60  90 

Qi  Qi 


Figure  3.  Diffusion  coefficients  Dj  and  D2,  in  s“*,  for  the 
variables  Qi  and  Q2  determined  by  the  curves  of  Figure  2. 
Qi  is  chosen  to  be  OO)  and  the  constant-Q2  curves  are 
labeled  by  the  value  of  logio(£'/0.2  MeV)  at  oio  =  90°. 


models,  the  wave  power  spectral  density  was  taken  to  be  the 
product  of  truncated  Gaussians  in  frequency  and  wave 
normal  angle  (tan  0)  in  each  of  three  different  local  time 
sectors.  In  the  nighttime  sector  model  (2300-0600  MLT), 
the  waves  lie  between  0  and  15°  latitude;  the  corresponding 
range  of  fpjfce  is  about  3.4  to  2.5.  The  waves  are  only 
present  between  15°  and  35°  latitude  in  the  prenoon  sector 
model  (0600-1200  MLT)  and  10°  to  35°  latitude  in  the 
afternoon  sector  (1200-1500  MLT);  the  corresponding 
ranges  of fpjfce  are  ~3.0  to  0.9  and  5.9  to  1.4,  respectively. 
The  tan  0  distribution  was  further  restricted  to  less  than 
0.9  times  the  resonance  cone  value,  and  resonance  harmonics 
up  to  n  =  ±5  were  included.  The  appropriately  weighted, 
combined  diffusion  coefficients  are  shown  in  Figure  1, 
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Figure  5.  Electron  flux  as  a  function  of  Oo  at  iE  =  0.5  MeV 
(left)  and  E  =  2  MeV  (right),  at  t  =  0  (blue)  to  f  =  1  day 
(red),  at  intervals  of  0.1  day,  corresponding  to  Figure  4. 


and  reflect  the  artificially  sharp  cutoffs  used  to  model  the 
frequency,  wave  normal  angle,  latitude,  and  local  time 
distributions  of  wave  power.  Note  that  the  cross  diffusion 
coefficient,  Da^p  ~  (iioio/^p/At),  can  have  either  sign. 

[10]  The  diffusion  coefficients  were  used  to  numerically 
trace  curves  of  constant  Q2  in  the  (oo,  E)  plane,  shown  in 
Figure  2.  (Since  computation  of  the  diffusion  coefficients  is 
quite  demanding,  equation  (8)  was  integrated  by  interpolat¬ 
ing  a  table  of  precomputed  values;  this  also  tended  to 
smooth  the  discontinuities.)  Diffiision  in  Qi(=  oq)  proceeds 
along  curves  of  constant  Q2,  which  are  roughly  curves  of 
constant  energy  (corresponding  to  the  approximation  of 
pure  pitch  angle  diffusion).  Diffusion  in  Q2  proceeds  along 
curves  of  constant  Qi,  which  reach  large  values  of  energy. 

[11]  So  far,  the  actual  values  of  Q2  have  not  been 
specified;  they  may  be  assigned  in  any  smooth  way  that 
labels  the  curves  on  which  they  are  constant.  From  Figure  2, 
it  is  natural  to  associate  with  each  constant-22  curve  the 
value  of  E  (more  precisely,  logio  [£'/0.2  MeV])  near  oto  = 
90°.  Intersections  of  the  constant-gi  curves  and  constant-22 
curves  establish  the  correspondence  between  coordinates 
(Qu  Q2)  and  (qo,  E).  To  evaluate  Di  and  D2  it  is  also 
necessary  to  evaluate  the  partial  derivatives  relating  the  old 
and  new  variables.  A  general  procedure  is  to  trace  curves 
with  constant  but  slightly  differing  values  of  Q2,  which 


E  (MeV)  E  (MeV) 


Figure  4.  Phase  space  density  0cft)X2n  Qi)  and  (right) 
foo,  £)  at  r  =  0,  0.1,  0.5,  and  1  day,  from  a  numerical 
solution  to  the  (2i,  Qi)  diffusion  equation.  The  lower 
boundaries  in  both  columns  correspond  to  £  =  0.2  MeV. 
The  initial  and  boimdary  conditions  are  discussed  in  the 
text. 


Figure  6.  Electron  flux  j-p^  fas  a  function  of  E  at  (left) 
Oo  =  30°  and  (right)  oo  =  70°  at  f  =  0  (blue)  to  t  =  1  day 
(red),  at  intervals  of  0.1  day,  corresponding  to  Figure  4 
(solid  curves).  Also  shown  are  results  at  f  =  0, 0.1,  and  1  day 
obtained  by  neglecting  cross  diffusion  Do^p  (dashed 
curves). 
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intersect  a  constant-gi  curve  at  slightly  different  values 
of  (Ji,  J^i).  This  yields  finite  difference  approximations  for 
dJ\ldQ2  and  dJ^ldQ^:,  and  similarly  for  derivatives  with 
respect  to  Qi.  These  determine  the  value  of  the  Jacobian  G 
along  with  dQ\ldJi  =  {dJ2ldQ^IG,  etc.  However,  with 
2i  =  Oo.  it  is  only  necessary  to  numerically  evaluate 
dEldQ2.  The  resulting  values  of  Di  and  D2  are  shown 
in  Figure  3. 

4.  Results  and  Discussion 

[12]  Using  the  values  of  Di  and  D2,  equation  (7)  was 
solved  numerically  in  the  (Qi,  Q2)  plane,  using  fully  explicit 
finite  differencing  for  simplicity.  To  qualitatively  model  the 
depleted  energetic  electron  levels  that  typically  follow  the 
main  phase  of  a  magnetic  storm,  the  initial  value  of 
differential  electron  flux  was  taken  to  be  j  =  exp  [— (JF  — 
0.2)/0.1]  sin  oio,  with  E  in  MeV,  and / = jlp^.  These  values 
were  held  fixed  at  the  lower  boundary,  namely  the  curve 
where  E{Q\,  Q2)  =  0.2  MeV,  while /=  0  was  imposed  at  the 
upper  boundary.  Zero  flux  was  also  enforced  at  the  loss 
cone  value  of  Qi  =  Oo,  and  dfdQ\  =  0  was  taken  as  the 
boundary  condition  at  Qi  =  90°. 

[13]  Results  after  1  day  are  shown  in  Figure  4.  Diffusion 
tends  to  transport  phase  space  density  to  higher  energy,  but 
the  outcome  of  the  competition  between  this  and  difl^ion 
into  the  loss  cone  could  not  be  anticipated  in  detail.  Pitch 
angle  profiles,  at  constant  energy,  are  shown  in  Figure  5.  At 
0.5  MeV  the  curves  evolve  roughly  as  expected  on  the  basis 
of  pitch  angle  diffusion  alone,  becoming  peaked  at  oo  =  90° 
and  falling  rapidly  near  the  edge  of  the  loss  cone.  However, 
at  2  MeV  the  profiles  develop  a  persistent  peak  around  35°, 
for  which  the  energy  diffusion  rate  is  largest.  Figure  6 
shows  one-dimensional  flux  vs.  energy  profiles  developing 
in  time.  Significant  levels  of  flux  develop  at  1  MeV  in  less 
than  half  a  day.  Energization  is  faster  at  cvo  =  30°  than  at 
70°,  as  expected  from  Figure  1,  and  reaches  ~2  MeV  during 
the  1  day  of  the  simulation. 

[14]  To  gauge  the  effect  of  cross  diffusion,  the  entire 
procedure  was  repeated  with  set  to  0.  The  results  are 
shown  in  Figure  6  as  dashed  curves  at  t  =  0,  0.1,  and  1  day. 
The  qualitative  behavior  is  similar,  but  for  small  Oo  neglect¬ 
ing  cross  diffiision  evidently  leads  to  an  overestimate  of 
energy  diffusion  (contrary  to  a  conjecture  by  Home  et  al. 
[2005]).  It  seems  likely  that  D^p/{Da^aJ)pp)  is  a  controlling 
parameter,  based  on  the  Green  function  expressions. 

[15]  Several  additional  physical  processes  must  still  be 
included.  Scattering  by  electromagnetic  ion  cyclotron 
waves  is  likely  a  strong  loss  mechanism  for  electrons 
above  1  MeV,  especially  combined  with  whistler  mode  hiss 
[Summers  et  al,  1998;  Summers  and  Thome,  2003;  Albert, 
2003].  Radial  diffusion,  particularly  when  enhanced  during 
storm  times,  is  a  potentially  crucial  source  of  energization 
[O’Brien  et  al,  2003]  and  also  serves  to  redistribute  locally 
accelerated  particles  [Green  and  Kivelson,  2004].  The 
framework  presented  here  should  help  make  it  possible 
to  realistically  combine  and  evaluate  Ae  effects  of  diffu¬ 
sion  in  all  three  adiabatic  invariants. 
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